Gravity v4 - Draft#1235
Conversation
…l errors in test values that appeared in test_flight_path_ode.py due to the use of altitude_rate.py and its english 32.2ft/s**2 to metric 9.81m/s**2 conversion that is now enforced by openmdao to be a conversion from 32.2ft/s**2 to 9.81456m/s**2
| 'load_factor': [0.25086942, 0.25086942], | ||
| Dynamic.Mission.ALTITUDE_RATE: [0.0, 0.0], | ||
| Dynamic.Mission.ALTITUDE_RATE_MAX: [-0.01815763, -0.01815763], | ||
| Dynamic.Mission.ALTITUDE_RATE_MAX: [-0.018143, -0.018143], |
There was a problem hiding this comment.
Why do these values change?
| # from gaspmain.f | ||
| # reli = reynolds number per foot | ||
| RELI = np.sqrt(1.4 * GRAV_ENGLISH_GASP * 53.32) * EM * np.sqrt(T0) / XKV # dynamic | ||
| RELI = np.sqrt(1.4 * grav_english * 53.32) * EM * np.sqrt(T0) / XKV # dynamic |
There was a problem hiding this comment.
This might be outside of the scope of this PR, but I think this equation looks like it is assuming air, and calculating the speed of sound and multiplying by mach number to obtain the V in reynolds number.
1.4: The ratio of specific heats for air.
53.32: The specific gas constant for air in English units.
Perhaps these values should be obtained from the atmopshere component so that it correctly calculates for gasses other than air.
| CM1 = 0.99850 # Center of mass (Earth)? Unknown | ||
| OC2 = 26.76566e-10 # Unknown |
There was a problem hiding this comment.
| CM1 = 0.99850 # Center of mass (Earth)? Unknown | |
| OC2 = 26.76566e-10 # Unknown | |
| CM1 = 0.99850 # Oblateness/Gravity exponent (accounts for earth not being perfect sphere) | |
| OC2 = 26.76566e-10 # Accounts for centrifugal acceleration due to Earth's rotation |
From conversation with ChatGSFC:
Constants Derivation for Geopotential to Geometric Altitude Conversion
CM1 (Oblateness/Gravity Gradient Term): 0.99850
- Represents the modified inverse-square law exponent for Earth's gravity.
- A perfect point-mass sphere yields an exponent of 2.0.
- Accounting for Earth's J2 harmonic (equatorial bulge/oblateness), the effective gravity
gradient exponent at standard latitude drops to ~1.9985. - Since the equation uses (R_earth / R)**(CM1 + 1.0), CM1 is defined as 1.9985 - 1.0 = 0.9985.
OC2 (Centrifugal Relief Term): 26.76566e-10 (approx. 2.677e-9)
- Represents the centrifugal acceleration due to Earth's rotation, which counteracts gravity.
- Derived from Earth's angular velocity (ω ≈ 7.2921e-5 rad/s) at the standard atmospheric
reference latitude (φ = 45° 32' 33"). - It is proportional to the centrifugal potential: ω² * cos²(φ).
(7.2921e-5)² * cos²(45.5425°) ≈ 2.676e-9 .
|
|
||
| CM1 = 0.99850 # Center of mass (Earth)? Unknown | ||
| OC2 = 26.76566e-10 # Unknown | ||
| GNS = 9.8236930 # grav_accel_at_surface_earth? |
There was a problem hiding this comment.
Since this function used a different constant before (we have now updated it to use the single gravity constant = 9.80665) then this function will evaluate differently. I think we are now accounting for the centrifugal relief twice. GNS is the old "newtonian gravity", which would be the gravity we feel if the earth isn't spinning. Our new grav_metric value already includes the effect of spinning and so we are now double counting it.
From discussion with ChatGSFC:
Why the Old Algorithm Works
In the old algorithm:
GN = GNS * (radius_earth / R) ** (CM1 + 1.0)
H = (R * GN * ... - Z * (R - Z / 2.0) * OC2) / g
It starts with GNS (9.8236930), which is the pure Newtonian gravity at the Earth's surface, ignoring rotation.
It decays that pure gravity with altitude to find the local pure gravity (GN).
It calculates the gravitational potential energy.
It subtracts the centrifugal potential energy explicitly with the OC2 term.
Finally, it divides by g (which is standard effective gravity, 9.80665) to convert the total potential energy into Geopotential Altitude (
Note: This exactly mirrors the official equations from the U.S. Standard Atmosphere 1976.
Why the New Algorithm is Bugged
In the new algorithm:
GN = grav_metric * (radius_earth / R) ** (CM1 + 1.0)
H = (R * GN * ... - Z * (R - Z / 2.0) * OC2) / grav_metric
It starts with grav_metric (9.80665). This number already has the surface centrifugal force subtracted from it.
It calculates the potential energy based on this artificially lowered gravity.
Then, it subtracts the OC2 centrifugal term AGAIN.
By doing this, the new algorithm penalizes the aircraft for Earth's rotation twice. As you go higher in altitude, this error will compound, resulting in an incorrect Geopotential Altitude.
There was a problem hiding this comment.
Maybe we should add a test for this utility and compare to a data table from somewhere to make sure this analysis is correct?
| R = radius_earth + Z | ||
| GN = GNS * (radius_earth / R) ** (CM1 + 1.0) | ||
| H = (R * GN * ((R / radius_earth) ** CM1 - 1.0) / CM1 - Z * (R - Z / 2.0) * OC2) / g | ||
| GN = grav_metric * (radius_earth / R) ** (CM1 + 1.0) |
There was a problem hiding this comment.
| GN = grav_metric * (radius_earth / R) ** (CM1 + 1.0) | |
| GN = (grav_metric + (OC2 * radius_earth)) * (radius_earth / R) ** (CM1 + 1.0) |
I believe this would correct the calculation, we might prefer to calculate the constant outside the loop for very marginal performance gains:
GNS = grav_metric + (OC2 * radius_earth)
2 participants
Summary
All gravity constants are now uniform and encapsulated in constants.py . The differences between FLOPS and GASP gravities have been eliminated.
Follow on to #1219, #1227, #1232.
Removes GRAV_ENGLISH_GASP and replaces with Mission.GRAVITY.
One reference to GNS as gravity in a standard python function was replaced with constants.GRAV_EARTH.
Related Issues
Partially resolves #1192
Backwards incompatibilities
None
AI Usage
None.